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Abstract 


Treatment of convective and pressure fluxes in the Euler and Navier-Stokes equations using 
splitting formulas for convective velocity and pressure is investigated. Two schemes — Controlled 
Variation Scheme (CVS) and Advection Upstream Splitting Method (AUSM) — are explored for 
their accuracy in resolving sharp gradients in flows involving moving or reflecting shock waves as 
well as a one-dimensional combusting flow with a strong heat release source term. For 
two-dimensional compressible flow computations, these two schemes are implemented in one of 
the pressure-based algorithms, whose very basis is the separate treatment of convective and pressure 
fluxes. For the convective fluxes in the momentum equations as well as the estimation of mass fluxes 
in the pressure correction equation (which is derived from the momentum and continuity equations) 
of the present algorithm, both first- and second order (with minmod limiter) flux estimations are 
employed. Some issues resulting from the conventional use in pressure-based methods of a 
staggered grid, for the location of velocity components and pressure, are also addressed. Using the 
second-order fluxes, both CVS and AUSM type schemes exhibit sharp resolution. Overall, the 
combination of upwinding and splitting for the convective and pressure fluxes separately exhibits 
robust performance for a variety of flows and is particularly amenable for adoption in 
pressure-based methods. 


1 . Introduction 

For the hyperbolic system of Euler equations, a number of high resolution schemes have 
been proposed to date. Most of these schemes are designed to satisfy the TVD property for scalar 
conservation laws and systems of equations with constant coefficients, whereby spurious 
oscillations are suppressed. Several different approaches can be found in literature such as the 
modified flux approach (scalar diffusion) of Harten [1], flux vector splittings [2, 3], flux difference 
splittings [4, 5], etc. All these schemes developed for the Euler equations can be directly extended to 
the Navier-Stokes equations. 



The main motivation behind the various approaches is to achieve high accuracy and 
efficiency in numerical computations, especially for complex flows that may involve strong 
convective effects, sharp gradients, recirculation, chemical reactions and / or turbulence models. 
Different schemes have different accuracy and efficiency characteristics. For example, flux vector 
splitting schemes are quite efficient and relatively simple, but produce excessive smearing. 
Moreover, the Steger- Warming splitting [2] produces glitches at points where eigenvalues change 
sign, such as sonic points. Van Leer splitting [3] is designed to remedy this, but it suffers from 
excessive numerical diffusion in viscous regions. Subsequent efforts have been made to reduce this 
diffusion (Hanel & Schwane [6]). On the other hand, flux difference splittings such as Roe and Osher 
splittings, have substantially lesser numerical diffusion. However, they too are known to yield 
inaccurate results in some simple flows. For example. Roe splitting produces non-physical 
“carbuncle” shocks in supersonic flows over blunt bodies [7]. 

Taking all the above factors into account, there is continued interest and ongoing effort in the 
development of new schemes which are robust in terms of accuracy as well as efficiency. Towards 
this end, one promising approach that is currently under investigation is the treatment of convective 
and pressure fluxes as two separate entities. Employing this idea, Thakur and Shyy [8-11] have 
developed a Controlled Variation Scheme (CVS) in which the convective flux is estimated using 
Harten’s second-order TVD scheme (modified flux approach) where the local characteristic speeds 
of the different equations are coordinated by assigning them the values of the local convective 
speeds; the pressure terms are treated as source terms and are central differenced or treated in a 
special manner by employing Strang’s time-splitting technique. Liou and Steffen [7] have also 
proposed a scheme, called the Advection Upstream Splitting Method (AUSM) which treats the 
convective terms and the pressure terms separately. In the AUSM scheme, the interface convective 
velocity is obtained by an appropriate splitting and the convected variable is upwinded based on the 
sign of the interface convective velocity. The pressure terms are also handled using an appropriate 
splitting formula. For both the CVS and AUSM schemes, the operation count is substantially smaller 
compared to Roe and Osher schemes since both CVS and AUSM do not entail the evaluation of 
either Jacobian matrices or intermediate states. 

In the controlled variation scheme (CVS) presented here, guided by the eigenvalues of the 
total flux as well as the individual convective and pressure fluxes, the treatment is as follows. The 
convective flux is fully upwinded whereas the pressure flux is split yielding contributions from 
upstream and downstream neighbors. Two different formulations which lead to different pressure 
fluxes, are investigated. The eigenvalues of the respective pressure fluxes are used to interpret the 
physical significance of the two formulations. It is shown that the formulation which is consistent 
with the physical mechanism that the convective fluxes get transported at the mean convection speed 
and the pressure signals propagate both upstream and downstream in subsonic flows, is perhaps the 
more desirable one. Two one-dimensional test cases — the standard shock tube problem and a 
longitudinal combustion instability problem previously investigated by Shyy et al. [12] — are used 
to demonstrate that the CVS and the AUSM scheme yield accuracy comparable to the Roe scheme. 
The results for the combustion instability problem, in particular, illustrate that the approach of 
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treating convective and pressure fluxes separately can indeed coordinate signal propagation even in 
the presence of source terms such as heat release. 

As far as the solution of the multi-dimensional Navier-Stokes equations is concerned, two 
broad categories of algorithms are in common use, namely the density-based and pressure-based 
algorithms. The attention of the present work is on the pressure-based methods. These types of 
algorithms, though originally developed for incompressible flows [13], can be easily extended for 
compressible flows [14]. The pressure-based algorithms treat the convective and pressure fluxes as 
two separate entities. However, to date, no theoretical foundation has been laid for incorporating the 
modem shock capturing schemes into these pressure-based algorithms for compressible flows. In 
the present work, the controlled variation scheme (CVS) developed earlier in the context of 
pressure-based methods for incompressible flows [10,11] is extended for compressible flows. The 
present work investigates the applicability of the recent pressure splitting formulas proposed by Liou 
& Steffen [7] in the context of their AUSM scheme into the pressure-based methods. It is 
demonstrated that via the CVS and AUSM type schemes, the pressure-based algorithms can indeed 
yield accurate simulations of complex compressible flows including high resolution of shock waves. 

2. Estimation of Fluxes for the CVS and AUSM Schemes 


The one-dimensional system of conservation laws for an ideal gas is the following: 

dW . dF _ ft 

dt d x 


where W = 


r Q ' 
m 

F = 

QU 

mu + p 


' Q u ' 
mu + p 

E 

< > 


(E + p)u 


Hu 


(la) 


(lb) 


Here, E is the total energy, E=g(e + u 2 /2), and H=E+p is the total enthalpy. A numerical scheme for 
Eq. (la) can be written as, for example 

w? +1 + - f^! /2 ) = W! - A(1 - B) - F”_ 1/2 ) (2) 

where X=At/Ax, F i± x / 2 are the numerical fluxes at the control volume interfaces, the superscripts n 

and n+1 represent time levels and 9 is a measure of implicitness of the scheme. We obtain explicit, 
fully implicit and Crank-Nicolson schemes for 9 - 0, 1 and 1/2, respectively. 

A recent approach is to treat convection and acoustic wave propagation as physically distinct 
(but coupled) mechanisms. The breakup of the total flux into convective and pressure fluxes can be 
done in at least two different ways as presented next. 


(a) Formulation 1: Based on Total Enthalpy 

One way of breaking up the total flux into convective and pressure fluxes is to treat the total 
energy flux (Hu) as part of the convective flux. Thus, the pressure flux consists of just the p term in 
the momentum flux: 
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'gu' 


'O' 

F = F c + FP = 

mu 

Hu 

^ J 

+ 

P 

0 



" ga " 


O' 

= M<P + F p = M 

gua 

gHa 

+ 

p 

0 


(3a) 

(3b) 


Such a breakup of the total flux has been used, for example, in the AUSM scheme of Liou and Steffen 

m. 


(b) Formulation 2: Based on Total Energy 


Another way of breaking up the total flux into convective and pressure fluxes is to treat the 
energy flux (Em) as part of the convective flux. Thus, the pressure flux now consists of p and pu 
terms: 



r QU' 


'O' 

F = F c + F p = 

mu 

Eu 

s a 

+ 

p 

pu 

W J 


(4) 


We first present the treatment of the convective flux for either of the above two formulations. 


2.1 Convective Fluxes 


2.1.1 Controlled Variation Scheme (CVS) 


As developed previously, the CVS utilizes the form of TVD type schemes — originally 
defined by Harten [1] — while defining the characteristic speeds in a different way [10]. The 
numerical convective flux F c . + using a first-order TVD scheme [1] can be written as: 


*7+1/2 


where Q is the convective dissipation function given by 

fc 2 


Qi +1/2 ” 0(^1+ 1/2) 


Ji( 


6 

\b\ , 


+ <5 


and 


A; 


i+1/2 


w = w i+l - w t . 


, + l/2 W } 

(5) 

if \b\ < 6 
if \b\ > <5 

(6a) 


(6b) 


The parameter <5 in Eq. (6) is used to eliminate the violation of the entropy condition for 
characteristic speeds close to zero [1] and b [ + ^ ^ is the local characteristic speed on the right 

interface of the control volume. 

Let w represent the dependent variable of each of the scalar conservation laws comprising 
the system (1) and let /represent the respective convective flux: 

/i+1/2 = fi + fi+ 1 _ ®( b i+1/2) ^i+l/2 W } ( 7 > 
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The local characteristic speed b 


is defined as follows: 


and 


^ i + 1/2 

A i+ 1/2 W 


* + 1/2 
Jj+l ~ /> 

, ^f+l/2 W 

V 

du 

k. 


if J m/2 W 5,4 0 


if J i+ i /2 w = 0 


(8a) 

(8b) 


In the present work, we employ the explicit scheme (0 = 0) for one-dimensional unsteady 
flow problems and the fully implicit scheme (i.e., 0=1) for multi-dimensional steady flow cases 
as the basis for development of the controlled variation scheme (CVS). For the latter, the implicit 
and highly nonlinear equations would require iterations at every timestep if a time-stepping 
approach to steady state is employed. If an infinite timestep is chosen to solve for steady state, as 
in the present study, the number of iterations required to achieve convergence will be very large. 
Consequently, some linearized versions of implicit TVD schemes have been devised [15-17]. We 
base our scheme on the linearized non-conservative implicit (LNI) scheme described by Yee 
[15,16], following which f i+l / 2 ~ /,_ 1 / 2 can be written as 



The superscripts « and n+1 signify the previous and current iteration levels at steady state, 
respectively. The above nonlinear equation can be linearized by dropping the superscripts of the 
coefficients of A f ± x / 2 w n + 1 from n+1 to n. This form can be shown to be TVD [16]. This form of 

the implicit scheme cannot be expressed in conservation form and thus it is non-conservative except 
at steady state where it has been shown that it does reduce to a conservative form [16]. 

In the CVS, the local characteristic speed b i+ x j 2 for the system (1) is defined as the local 

convective speed: 

^x+l/2 = \( u i + u i+ 1) < 10 ) 

Also, in the CVS, the parameter <5 is used to regulate the amount of numerical dissipation in the 
scheme. 


2.1.2 The AUSM Scheme 

We now briefly present the treatment of convective flux in the AUSM scheme proposed by 
Liou and Steffen [7] which also treats the convective terms and pressure terms separately. The 
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( 11 ) 


numerical convective flux, for the AUSM scheme, at an interface is written as: 

F i + 1 /2 = M i + 1 /2 (*i + 1 / 2 )^ 
where M i+ x j 2 is the interface convective velocity and <P is the convected variable. 

The interface convective Mach number is expressed as the sum of the split values of the 
positive and negative contributions from the left and right states of the interface: 

M i + 1/2 = (^1+1/2)^ + (^i+1/2) 


' R 


( 12 ) 


For first-order accuracy, the left (L) and the right (R) states on the (i+1/2) interface are obtained by a 
first— order extrapolation of the nodal values of the left and the right neighbors of the interface: 

K + , KK =M ' +i (13) 

The spitting chosen here (based on the van Leer splitting for the Euler equations) is the following: 

if IMI < 1, 

(14) 

otherwise 


M ± = 4 


± ±(M ± l) 2 , 

\{M ± IMI), 


The convected variable is upwinded, depending on the sign of the interface velocity, as 


follows: 


K+l /2 )l 

i l \ 

1-^ 

+ 

N> 

IV 

p 



(15) 

(^•+1/2)^ 

otherwise 



(^ t+ l ^L/R 

For first-order accuracy, a first-order extrapolation using one upwind nodal value is employed: 
(^1+1/2)^ = (®i+l/ 2 ) R = ^«+i ( 16) 


2.1.3 Formulation of the Fluxes of the CVS and AUSM Schemes 


Wada and Liou [18] have proposed a version of the AUSM scheme based on flux difference 
splitting, labelled as AUSMD. It is interesting to note that the estimation of the convective fluxes 
in the momentum equations by the AUSMD scheme is identical to that of the linearized CVS scheme 
for steady state computations, except the interface convective velocities are estimated differently 
in the two schemes. The net convective flux along the x-direction, for example, using the CVS 
scheme, Eq, (9), with <5=0, can be expressed as follows: 




CVS 


2^i+l/2( M i+l/2 
2^<-l/2(“ u i- 1/2 




(“«+] 


£i+I/ 2 l M i+l /2 




U i—\/2 “ 1 - 1 / 2 ) ( U i U i~ l) 


(17) 
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In the AUSMD scheme, the convective flux is expressed as follows: 

1 AUSMD 

M.+ 1/2 ” iQu2) i- 1/2 = 2l (eU) i+ld U i + “i+ 1) " M.+ l/J (“«+l “ “»)] 

~ ^ (Q“)i-i /2 ( u i - 1 + M i) " |(e«) f -i/ 2 | (“« “ u «-l)] 

= ^i+l/ 2 (“i+l/2 - K+ 1 /J) K+l - M <) - 2 ^ ,— 1 / 2 ( _ 1/2 " K-I/2I) (“i - M «-l) 

+ M < (^i+l/ 2 w j+l /2 ~ l/ 2 W i- I/ 2 ) (18) 


From Eqs. (17) and (18), it can be seen that 

-.AUSMD 

e“ J ) i+1/2 - fe“Vi/2j = [(e« 2 ) i+1/z - (<?“V,/ 2 

+ U i (^i+ ]/ 2 U l+ 1/2 1 / 2 M j — I/ 2 ) (19) 

The last term in the above expression is nothing but the nodal value of the dependent variable 
multiplied by the net mass flux in the x-direction. A similar expression results from the convective 
fluxes along the y-direction. Thus the difference between the numerical convective fluxes between 
the CVS and the AUSMD schemes is the net mass flux term integrated over a control volume, which 
must be zero at steady state (from the continuity equation). Thus, for steady state applications, the 
two flux estimations are identical. The only difference is the method of estimation of the interfacial 
velocities — the CVS scheme just averages the nodal point values (Eq. 10) whereas the AUSMD 
scheme uses splitting based on the local Mach number (Eq. 12). 

2.2. Treatment of Pressure Flux 

As mentioned earlier, the pressure flux is treated separately in both CVS and AUSM type 
schemes. Different approaches can be taken as described next. 

2.2.1 Splitting of the p term 

The p term in the momentum equation for both the formulations given by Eqs. (3) and (4) can 
be treated by splitting as follows: 



For first-order accuracy, the left (L) and the right (R) states on the (i+1/2) interface are obtained by a 
first-order extrapolation of the nodal p values of, respectively, the left and the right neighbors of the 
interface (Fig. 1). 

The splitting of pressure can be achieved in a manner similar to the van Leer splitting for the 
fluxes of the Euler equations [3], The van Leer splitting is based on the requirements that the split 
fluxes as well as their first-derivatives be continuous and that the split fluxes be polynomials of the 
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lowest possible degree. This leads to a splitting of the fluxes in terms of factors (M ± 1 ) . Liou et 
al. [19] have suggested a similar splitting for pressure: 

p = a[(M + l) 2 - (M — l) 2 ] + P[(M + l) 2 + (M - l) 2 ] (21) 

= 4 aM + 2 P{M 2 + 1 ) 

By choosing 2/S=p (which is true for M= 0), we obtain 

a = — ^pM (22) 

and hence p can be written as 

p = l(M + 1) 2 (- M + 2) + |(Af + 1) 2 (M + 2) (23) 

Thus, as suggested by Liou and Steffen [7] and Liou et al. [19], the following splitting is employed: 

j(M ± 1) 2 (2 =F M), 

± ^ 

P ~ * |(M ± 1M1)/M, 


iflMI < 1, 
otherwise 


(24) 


Thus, for supersonic flow, the above formulation leads to full upwinding of pressure, namely, 

p. + 1/2 = Pi , if IMI > 1 and M > 0 (25) 


The splitting takes place only for subsonic flow where contributions from both upwind and 
downwind neighbors are taken into account: 


Pi+ 1/2 



if IMI < 1 


= i (M, + lf(2 - M,) Pi + i(M j+1 - l) 2 (2 + M j+1 ) Pi+1 (26) 


This is consistent with the fact that pressure signal propagates only in the upwind direction for 
supersonic flows and in both upwind and downwind directions for subsonic flows. The pressure flux 
at the interface can also be interpreted as the following: 



P( + 1/2 + 2^i*+l/2 

(27a) 

with 

= ./ 2 (3 - "? + ./a) 

(27b) 


where Q p ^ is the numerical diffusion introduced by the pressure splitting to the central difference 
flux and M . + j j 2 is the interface Mach number. It can be observed that as M . + j j 2 increases from 0 to 
l t QP+ j ^varies from 0 to 1 smoothly, changing the nature of the pressure splitting from central 
differencing (for M i + 0) to full upwinding (for M. + = 1) tn a continuous manner. 
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Eigenvalues 


We next present a brief analysis of the eigenvalues associated with the convective and 
pressure fluxes. It must be stated at the outset that one has to examine the eigenvalues of the total 
combined flux (convective and pressure) in order to interpret the true nature of signal propagation in 
the gas dynamic system. However, the following analysis gives an idea of the nature of the 
convective and pressure fluxes in the two formulations given by Eqs. (3) and (4). In particular, from 
the viewpoint of operator splitting, the following analysis will be relevant to the individual 
components of convective and pressure fluxes, respectively. Also, this breakup of the total flux is 
expected to be more critical for the CVS since in this scheme, the coupling between the convective 
and pressure fluxes is not explicitly coordinated as a function of the local Mach number, as in the 
AUSM scheme. For the convective and pressure fluxes given by Eq. (3), the Jacobians of the 
convective and pressure fluxes are given by the following 

1 

2 u 

y § - §(y - i )« 2 


A c = M£ = 
dW 


0 


—u 


-r fu + (y 


- l)ir 


0 

0 

yu 


(28a) 


A p s 


dFP _ 
dW 


0 

u 2 

(y - D*- 

o 


0 

- (y - 1)« 
o 


0 

(y - i) 
o 


(28b) 


where y is the ratio of specific heats of the gas. The eigenvalues of the above Jacobians can be found 
by solving the equations 

A c — A c I = 0 A p — A p 1 = 0 (29) 

where /is the identity matrix, resulting in the following eigenvalues for the convective and pressure 
fluxes: 


A c = 

'O' 
2 u 

A p = 

0 

- (y - i)« 


yu 

K * 


0 


which seem to indicate that the convective flux has an upwind character and that the pressure flux has 
a downwind character only. In accordance with the physical characteristics, it is desirable that the 
convective fluxes are completely upwinded and pressure fluxes are split based on the local Mach 
number. It is this thought which prompts us to investigate Formulation 2 (Eq. 4) which is perhaps 
more consistent with the numerical treatment of the individual (convective and pressure) fluxes. 

If the three equations — continuity, momentum and energy — are looked upon as a 
collection of three scalar conservation laws, as in the case of sequential solvers [8,9], then the 
characteristic speeds that one obtains for the convective terms in the three equations are 
( u, u, yu ). Likewise, for the pressure terms in the three equations, the characteristic speeds are 
( 0, -(y - 1 )u, 0 ). 
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2.2.2 Splitting of the pu term 

For formulation 2, given by Eq. (4), the pressure flux has the pu term in the energy equation, 
which can also be split in a manner similar to the p term in the momentum equation. As suggested by 

Liou et al. [19], the pu term can be expressed as consisting of (M ± 1) factors occurring in the 
splitting for p and a quadratic function in u: 

pu = ^(M + 1) 2 (a« 2 + 2 Bua + Ca 2 ) - ^(M - 1) 2 (Aw 2 - 2Bua + Co 2 ) (31) 


= (A + B)gu 2 + (B + Qgua 2 

in 

A + B = 0 
P 


The qu 2 term can be eliminated by enforcing the condition 
to obtain 


B + C = 


ga z 


(32a) 


(32b) 


Thus, a family of infinite choices for splitting pu are possible based on the parameter B. From the 
consideration of the total energy flux in the Euler equations, van Leer has proposed the following 
choice of B : 

h/a 2 


B - 


1 + Ihja 2 


where h = H/g. The simplest choice, as proposed by Hanel etal. [20] is 

5 = 0 


(33) 


(34) 


Thus, as suggested by Liou et al. [19], in the present study, the following splitting is employed: 

± |pa(M ± 1 ) 2 (Am 2 ± 2 Bua + Ca 2 ), if IM1 < 1, 

± \M\)/M, otherwise 


pu ± = < 


(35) 


Eigenvalues 

For the convective and pressure fluxes given by Eq. (4), the Jacobians of the convective and 
pressure fluxes are given by the following 

0 

,2 


A p = m = 
A ~ ew 


ac — dFl = 
A ~ aw 


o 

(y - i)f 


— U 

~Q U 


1 

2 u 


0 

0 


I “ 


0 


- (y - l)u 


(36a) 


(y 


-«(-f" + “ 3 ) <y ~ »(f - f“ 2 ) 


0 

(y - i) 
(y - i )u 


(36b) 


The eigenvalues of the above Jacobians are: 
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( 37 ) 


which indicate that the convective flux has an upwind character and that the pressure flux has both a 
downwind and an upwind character. The eigenvalues of the pressure flux in this formulation also 
suggest that the speed of propagation of pressure signals is dependent on the acoustic speed. Thus, it 
appears that Formulation 2 is more consistent with the dynamics of the system if one treats 
convection and acoustic wave propagation as two separate entities. Once again, it must be stated that 
Formulations 1 and 2 are expected to make more difference for the CVS than the AUSM scheme due 
to the reasons stated earlier. 

If the three equations in the Euler system (continuity, momentum and energy equations) are 
looked upon as a collection of three scalar conservation laws, then the characteristic speeds that one 
obtains for the convective terms in the three equations are ( u, u, u ). Likewise, for the pressure 
terms in the three equations, the characteristic speeds are ( 0, — (y - l)w, (y - l)w ). 



2.3 Roe Scheme 


In the present paper, the performance of the CVS and AUSM schemes is compared with that 
of the Roe scheme (which treats the convective and pressure fluxes as one combined flux and is 
considered as the most accurate of the approximate Riemann solvers) for some one-dimensional test 
cases. We briefly present the fluxes of the Roe scheme to contrast those of the CVS and AUSM. The 
Roe scheme [4, 21] is based on a characteristic decomposition of flux differences: 

F „ M - F, = A (*,+,„) K +1 - w,) <38) 

where A^W i+ 1/2 ) is the linearized Roe-averaged matrix . The numerical flux of the Roe scheme at 
an interface for each of the three equations is given by 

W = lb + fi* >) - ill U'W*’® (39) 


where A and R are respectively the eigenvalues and eigenvectors of a(w. + | /2 j: 


2, = u — c 


X 2 = u 


X 3 = u + c 


(40a) 


/?,=-% u-c 
2 c H — uc 


Ro = =2 


R, = u + c 


2c \H + uc 


(40b) 


Here the A I - +/ / 2 VV ’ S are the characteristic variables which represent the magnitudes of the waves at 
the interfaces and are given by: 
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(41a) 


where Au = A i+] j 2 u 
root of densities: 


zdw (1) = ^7 — Au 

QC 

Aw (2) = Aq (41b) 

c 

AwW = i?r + Au (41c) 


Qc 


u t , etc. All the variables are weighted by the ratio of the square 

^i+i/2 ” 

(42a) 

U i+ 1/2 = ^i+1/2 Qi 

(42b) 

= ^i+l/2 U i+l + U i 

U i+ 1/2 — J? , + 1 

A j+1/2 ^ 1 

(42c) 

= _ R i+l/2 H i+l + H i 

"i+ 1/2 — P , + 1 

A i+l/2 + 1 

(42d) 

“m/2 = O' - D^+i/z ~ 2“ 2 ) 

(42e) 


2.4 Numerical Dissipation of the Various Schemes 


The net interfacial flux for each of the schemes discussed in the previous sections can be 
expressed as follows: 


CVS: 

F i+ M2 =\[ F i + F Ul ~ Q i+ Ml ^ + 1/2 W ] + F i +1/2 

(43a) 

AUSM: 

F i +1/2 = 2 + 1/2( + ^f+l) _ |^i+ I/2I ^i+1/2^] + F i+ 1/2 

(43b) 

Roe: 

F i+\/2 = 2 [ + F i+ 1 — | -^i+1/2 | ^i+1/2^] 

(43c) 


In all of the above schemes, the last term in the square brackets represents the numerical dissipation 
added to the central difference scheme. A significant contrast among the above schemes is that the 

Roe scheme involves the computation of the linearized matrix j/ 2 )> unlike both CVS and 

AUSM schemes. In the AUSM scheme, the coupling between the numerical convective and pressure 
fluxes via the splitting formulas for both velocity and pressure at the control volume interfaces is 
expected to coordinate signal propagation, thus yielding no spurious oscillations. In the CVS, too, 
there is such a coupling but perhaps to a lesser degree since interfacial velocity is directly estimated 
by two-point averaging. However, the parameter 6 in the CVS serves to regulate the amount of 
numerical dissipation in order to suppress spurious oscillations should they occur. 
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3. Extension to Second-Order Spatial Accuracy 


The net flux for each of the schemes discussed in the previous sections is spatially first-order 
accurate. One can extend the net flux formally to second-order accuracy by employing a variable 
extrapolation or a flux extrapolation approach. The latter is chosen here. 

Let the total first-order flux at an interface be given by f i±l j 2 where the various quantities at 

the interfaces are defined by first-order extrapolations. We define f t = /(«,) as point-valued fluxes 
at cell centers. 

A general expression for a higher-order flux at the interface can be written as [22]: 



For x = 1, we obtain the second-order central difference scheme whereas x = -1 yields the fully 
upwind second-order scheme. Choosing x = -1, the second-order interface flux becomes: 


1/2 ~fi+ 1/2 + 2 fi-l/l) + (fi+1 ^+ 3 / 2 ) 


(45) 


which can be expressed with the use of limiters (for suppression of oscillations), as follows: 

f-l\, 2 = f i + 1/2 + 2 ^<-1/2) ■ {ft ~ fi - 1/2) - \ ^r + 3 / 2 ) • (/;•+! fi + 3/2) ( 46 ) 

where 


r 


+ 

i+l/2 


fi 


+ 2 


fi + ) 


fi+2H 
^i + 1/2 


r i+ 1/2 


fi-\ fi- 1/2 
fi ~fi+ 1/2 


(47) 


The function W(r) is the flux limiter mentioned above. The minmod flux limiter has been employed 
in the present study, namely: 

Minmod limiter : ^(r) = max[0, min(l, r)J (48) 


Other limiters such as van Leer’s monotonic limiter or Roe’s superbee limiter can also be used. 

For the implicit version of the CVS, the second-order interfacial flux can be written as 
follows: 

4 + 1/2 =f i+ 1/2 + 4 ^(^- 1 / 2 ) • [fi -fi- 1 + Qi-l/2 A i-l/2 w ) 

- | 4 r '"+ 3 / 2 ) • {fi+2 - fi+ 1 - Q i + 3/2 A i+ 3/2 w ) ( 49 ) 

which can be further simplified, using Eq. (8a), as follows: 

4”/2 = /i+ 1/2 + + Q <-y fl) ' ("'< " w i-i) 

+ ^(' 7 * 3/2) ■ j- *,+3/2 + 2, +3/2 ) • (w i+2 - W, + ,) (50) 

Similarly, can be written as 
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4-1/2 = f i- 1/2 + 4^( r «-3/2) • K-3/2 + <2i-3/ 2 } * ( W i~ 1 “ W i-2> 
+ ^r+1/2) * {- ^ + i /2 + e (+ i /2 } • K+i - ";) 


( 51 ) 


Finally, from Eqs. (7), (50) and (51), we get the following net flux for the linearized implicit 
version of the CVS 

•O 71 

\ n + 1 


p.) _ p.) _ I 


j + 1/2 J i- 1/2 


*>,+ 1/2 - Gi+1/2 1 + A^(^i+I/2) ’ h+i - W f )" 


+ i 


- *,- 1/2 - fif-,/2] [* + M r * + -l/2)]j ( W * ” W *-l)” + 1 

> n 

ty[ r i+ 3 / 2 ) [“ ^i+3/2 + 2, + 3 / 2 j ( w i +2 - w i+l) 

n 

x p[^ r i'- 3 / 2 ) [ ^i- 3/2 + ^- 3 / 2 ] ( w i-l — — 2 ) j 


(52) 


4. Results of One-Dimensional Test Cases 

Two one-dimensional cases are investigated in the present study. The first one is the standard 
shock tube problem. The second problem involves the simulation of a longitudinal combustion 
instability which involves thermoacoustic coupling due to the interplay between pressure 
oscillations and periodic heat release in the combustor. 

4.1 Shock Tube Problem 

The shock tube problem presented here has been previously investigated in Thakur & Shyy 
[8,9]. The initial conditions on the left and the right of the diaphragm are reported there. In the 
present work, we study the CVS and AUSM schemes in terms of their capacity to resolve 
discontinuities such as shock waves and contact surfaces. The total length of the tube is 14 units with 
the initial location of the diaphragm in the middle of the tube. 141 grid points are used and the value 
of X = At/Ax is 0. 1 , Results are presented in the form of total energy profiles after 200 time steps. For 
all the results presented in the following, the fluxes of all the schemes employed have been 
extrapolated to second-order using the minmod limiter. 

Figures 2 and 3 show the total energy profiles obtained with Formulations 1 and 2, 
respectively (as classified in Section 2), of the CVS, using two values of the parameter <5 which 
regulates the amount of numerical dissipation. It can be observed that Formulation 1, which treats 
only thep term as part of the pressure flux, yields a slight overshoot near the shock location, as seen in 
Fig. 2(a). For both the values of d. Formulation 2 yields solution profiles which are qualitatively 
better than those obtained with Formulation 1 , consistent with the interpretation of the eigenvalues 
of the Jacobian matrices for the two formulations, discussed in Section 2. Comparing these profiles 
with those obtained using the second-order Roe scheme, shown in Fig. 7, it can be observed that the 
CVS, especially with Formulation 2, yields accuracy comparable to the Roe scheme. 
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The solution profiles obtained with the second-order AUSM scheme are shown in Fig. 4. 
Again, both formulations for the convective and pressure fluxes are investigated. It is seen that the 
AUSM scheme yields results with Formulation 1 which are comparable with those obtained with 
Formulation 2. A possible explanation is that the AUSM scheme uses splittings for both the 
convective interface velocity and the pressure flux along with upwinding for the convected 
variables. Thus, the u-velocity of the eigenvalues in Section 2, for the AUSM scheme, already 
exhibits directional bias according to the local Mach number, which is not the case for the CVS. The 
AUSM scheme, like the CVS, also yields results of accuracy comparable to those obtained with the 
Roe scheme (Fig. 7). 


A Note on the Implementation of Pressure Splittings 


As far as the splitting for the pu term is concerned, numerous choices are available based on 
the parameter B, as discussed in Section 2.2.2. The splitting suggested by van Leer based on Eq. (33) 
and by Hanel et al. [20] based on Eq. (34) have both been implemented for the shock tube problem. 
They yield virtually the same results. 

The pressure splittings given by Eq. (24) for Formulation 1 and Eqs. (24) and (35) for 
Formulation 2 of the CVS can be implemented in different ways. One method is to use Eq. (26) for 
the p term where the nodal Mach numbers are used, and a similar equation for the pu term; this is the 
implementation employed for the results shown in Figs. 2 and 3. This is the original splitting for the 
pressure flux employed in the AUSM scheme of Liou & Steffen [7]. The other method is to use the 
interface Mach number, e.g., Eq. (27) for the p term. This expression for the pressure flux is similar 
to the one employed in the CUSP (Convective Upwind and Split Pressure) scheme of Jameson [25], 
The energy profiles obtained using this implementation for pressure splitting in Formulation 1 of the 
CVS are shown in Fig. 5. Although these two implementations appear to be very similar, they yield 
very different results, as seen from Fig. 2 and Fig. 5. The first implementation appears to be quite 
robust whereas the second implementation yields spurious oscillations in the entire solution profile. 
These oscillations can be reduced by increasing the amount of numerical dissipation (by increasing 
<5), but this is accompanied by an overall smearing of the solution profile, as seen in Fig. 5(b). 

The two above-mentioned implementations for the splitting of pressure flux have also been 
investigated for the AUSM scheme. Fig. 4 shows the results obtained with the first implementation 
based on the nodal Mach numbers which, as mentioned earlier, is the original pressure splitting 
employed in the AUSM scheme (Liou & Steffen [7]). Fig. 6 shows the results obtained with the 
second implementation of the pressure flux which is similar to the CUSP scheme (Jameson [25]). It 
can be seen that this approach is much more sensitive than the CVS and yields substantial oscillations 
when pressure splitting based on the interface Mach number is employed. Combining the various 
aspects discussed, including the choice of the dependent variables as well as the estimation of the 
variables on the control surfaces, it is clear that the detailed implementation is often as critical as the 
overall formulation of both the CVS and the AUSM schemes. 
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4.2 Longitudinal Combustion Instability Problem 

This test case has been devised by Shyy etal.[ 12] to investigate the interaction of convection 
and a source term in the form of heat release. It involves pre ssu re oscillations in a one-dimensional 
model of a combustor which are sustained by the oscillations of heat release. The heat release in the 
combustor is specified using a simple model based on some experimental observations. The details 
of the heat release model have been presented in Shyy et at. [12]. The interaction of the heat release 
source term and the convective and acoustic mechanisms in the system can lead to non-physical high 
frequency oscillations in some solution profiles. These numerical oscillations, besides being 
fundamentally undesirable from a numerical accuracy point of view, may lead to an instability of the 
computation or even trigger nonlinear instabilities in the system. In this regard, this problem is quite 
a stringent test case for any numerical scheme which seeks to coordinate convection and acoustic 
wave propagation with source term effects such as heat release effects. The value of X = At/Ax is 0.03 
and the simulation is carried on for 2x 10 6 time steps. Results are presented in the form of the 
following: 

(i) pressure and temperature mode shapes plotted at the last ten instants which are 10 4 
time steps apart, shown in Fig. 8, and 

(ii) pressure and heat release time series at the location x=0.75 for the last 4 x 10 4 time steps, 
shown in Fig. 9. 

The following schemes are used: 

(a) CVS with 6=0.0 and pressure splitting with Formulation 2 

(b) Second-order AUSM scheme 

(c) Second-order Roe scheme 

For the CVS, central differencing of the pressure flux along with a lower amount of damping (6=0.0) 
leads to spurious oscillations in some of the pressure mode shapes (see the results in Shyy et at. [12]). 
These oscillations can be suppressed with a sufficiently high dose of dissipation by increasing the 
value of 6 to 0.8 or higher [12]. This extra dissipation, of course, leads to a smearing of solution 
profiles which is manifested in the form of reduced magnitudes of all the mode shapes . This can be 
resolved by treating pressure as a source term along with heat release and imparting a special 
treatment such as Strang’s operator splitting method (Thakur & Shyy [8,9]). Such a special source 
term treatment yields an improved accuracy with no spurious oscillations [12]. However, this is at an 
increased computational expense which can be avoided by splitting the pressure flux in an 
appropriate manner as discussed earlier in this work. 

As seen from Figs. 8 and 9, the results obtained with the CVS along with splitting of the 
pressure flux, employing Formulation 2 of the convective and pressure fluxes are quite satisfactory 
and comparable in accuracy to those obtained with the second-order Roe scheme. The C V S is able to 
coordinate the interaction of pressure oscillations and the heat release in an appropriate manner and 
no spurious oscillations are observed. The CVS yields pressure mode shapes with slightly larger 
magnitudes compared to those resulting from the second-order AUSM and Roe schemes. 

The second-order AUSM scheme also yields results comparable in accuracy to those 
obtained with the second-order Roe scheme. However, slight oscillations are observed in mode 3 
(Fig. 8) of the pressure mode shapes in the region given by x=0.6 to x=0.8. Overall, the CVS, the 
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AUSM and the Roe schemes yield results very comparable in accuracy. However, one can observe 
from the mode shapes of the variables (especially pressure) that the results obtained by these three 
schemes are not exactly in phase. This is to be expected because there is a difference in the dispersive 
and dissipative characteristics among the schemes, which, however small, will result in a slight 
difference in the phase characteristics of the solution, especially after 2 X 10 6 time steps. 

5. Pressure-Based Solver for Multi-Dimensional Flows 

The pressure-based method employed for steady state solutions of the Navier-Stokes 
equations in the present study is in the spirit of the well known SIMPLE (Semi— Implicit Method for 
Pressure-Linked Equations) algorithm [13]. This algorithm was originally developed for 
incompressible flows and has been successfully extended to solve compressible flows by modifying 
the pressure correction equation to take the effect of variation of density on pressure into account 
[14, 23]. Additionally, the algorithm has been extended to body-fitted curvilinear coordinates in 
order to handle arbitrarily-shaped flow boundaries [14]. 

The present pressure— based algorithm employs a control volume approach and uses a 
staggered grid for the velocity components and the scalar variables like pressure, in order to avoid 
a checkerboard pressure distribution. The algorithm solves the governing equations in a sequential 
manner. The velocity components are computed from the respective momentum equations. The 
velocity and the pressure fields are corrected using a pressure correction equation which is derived 
by manipulating the continuity and the momentum equations. The correction procedure leads to a 
continuity— satisfying velocity field. The whole process is repeated until the desired convergence is 
reached. 

Both the CVS and AUSMD schemes are implemented in the present pressure-based 
algorithm. As pointed out in Section 2.1.3, the convective fluxes of the CVS and AUSMD schemes 
are identical for steady state applications except for the computation of the interface convective 
velocities. 

5.1 Momentum Equations 

For a two-dimensional fluid flow problem, the momentum equations in body-fitted 
coordinates can be written as follows: 

+ + = - «**») ] + + ?3 ^) ] + 

° r + ^ p )] (53) 

where x.,x v , and y v are the metrics associated with the curvilinear grid. A typical control volume 
is shown in Fig. 10. J (the Jacobian of the transformation), q v q 2 and q 3 are given by the following 

J = x^y v - yr,Xr, (54a) 

<7l = x y 4" y\ 
q*l ~ x £ x *i 4 " )*£ y*] 
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(54b) 

(54c) 



<?3 = x \ + y \ 


(54d) 


The contravariant components of velocity are given by: 

U - uy n - vx v 


(55a) 

(55b) 


V = VX| — uyg (55b) 

In the above, <p = u or v and fj. is the physical viscosity of the fluid. Integrating the above over the 
control volume and arbitrarily taking A£ — Arj — 1, we get 

[(p0)/> “ + [(.QU<t>) e ~ (QU<p) w } + [(Qv<f>) n - (pv^) 5 ] = 

^(<? 10 £ “ ^2^7)] - [7(^1^ ~~ ^2^7)] + [ 7 (” < 720 | + 93^7)] “ [^(” < 720 $ + ^3^7) 

+ {~yriP£ + y$Pri) or (* 7 P% ~ x ^Pv) ( 56 > 

where the superscript 0 represents the previous time level. For more details, see Shyy [14]. 

We now formulate the controlled variation scheme (CVS) by formally extending the scheme 
from one dimension to the present two-dimensional case. It should be recalled that the convective 
fluxes of the CVS and AUSMD schemes are identical for steady state applications except for the 
computation of the interface convective velocities. Using the form of the CVS presented in Eq. (52) 
independently along the x- and y-directions, Eq. (56) can be expressed, for steady flow, as follows: 

*u+i/J&e “ $p) ~ C U,i^p ~ 0w) + CJ +1/2 (<p N - <pp) - C"_j/ 2 (0/> “ 0s) 


— D n 
U i + 1/2 


{<Pe ~ 0J>) ~ D”- 1/2(0/* 0w) 


+ D j+\/2^N - 0p) “ 1/2(0^ + * 


where the various coefficients are given by 


± 1/2 — 2 1/2 ± ± 1/2 ^± 1 / 2 ] ^ + 2 ^( r «'± l /2 

c ,±i /2 = (j«i). ± 1/2 s (y«i) I ■ etc - 


(58a) 


(58b) 


and b i+ ,y 2 , etc. are the interface velocities. S is the source term consisting of pressure gradient and 

viscous cross— derivative terms as well as the remaining higher— order contributions from the 
convective fluxes. Note that the subscript i is used to denote the ^-direction and j to denote the 
^-direction. 

Using the conventional notation for the SIMPLE algorithm [13], Eq. (57) can be written as 
Ajjpp — A£0£ + + Af/p N + A<ff> s + S (59) 

where 
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[^'Ll/ 2 "2 e,+l/2 f 

[^ 3 L ./2 " 2 ®' +,/2 ■ 

Ap^ + A e + A 5 + Aft 


b i— 1/2 

^j+1/2 

_ ^y- 1/2 
b j+ 1/2 


- a ( - 1/2 ] [i +M^- 1/2 )] 

-< 2 1 + i / 2 ] [* +^ r + i / 2 )] 

- 0 ,_ 1 / 2 ][ i +^( o + - 1/2 )] 

-e ;+1/2 ] [l +^(rf +l/2 )] 


- 4 01 + 3/2 ' V'( r i+3/2) • { ^,- + 3/2 

+ 5 0 j — 3/2 • ^(^ + _ 3/2 ) ’ ( ^f-3/2 

- \ 0; + 3/2 ‘ W0"+3/2) • {- ^y+3/2 

+ 4 ’ 0/— 3/2 ' ^(^ 3 / 2 ) ' { *7-3/2 




1/2 


+ Gi+3/2} ‘ (0« + 2 ~ 0i+l) 
+ 2,- 3 / 2 } • (0,-1 - 0«-2) 
+ G7+3/2} • (0y+2 - < Pj+ 1) 
+ Qj - 3/2} * (0/-i _ 07-2) 




1/2 


(60a) 

(60b) 

(60c) 

(60d) 

(60e) 


(60f) 


+ P* 


In the above, P*represents terms involving pressure. It can be observed that the above form is 
spatially a five-point scheme along both directions which can be conveniently solved using the ADI 
method along with a tridiagonal matrix solver. Also, it should be noted that the coefficient matrix 
has a dominant diagonal. 

For the boundary control volumes, first-order numerical fluxes are employed to obtain the 
coefficients A E , A w , etc. 

5.2 Pressure Correction Equation 

The equation of continuity is represented indirectly by the pressure correction equation since 
pressure, not density, is the primary variable in the SIMPLE algorithm [13]. The pressure correction 
equation consists of mass flux terms ( qU )* and (pV) as well as correction terms on the 
pressure-correction control volume interfaces along the £- and rj — directions [14]. The superscript 
(*) indicates that these mass fluxes do not satisfy the continuity equation during the course of 
iteration (hence the need for pressure and velocity correction). The velocity correction terms are 
related to pressure corrections. Additionally, for compressible flows, the variation of density must 
also be accounted for via density corrections which are also related to pressure corrections through 
the equation of state. This changes the nature of the pressure correction equation from a pure 
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diffusion equation (for incompressible flows) to a convection-diffusion equation (for compressible 
flows). For details the reader is referred to Shyy [14, 23]. 

Conventionally, the estimation of the mass flux terms is done by utilizing the normal velocity 
components located at the faces of the pressure correction control volumes (due to the staggered 
grid) and by upwinding the density based on the direction of the interfacial normal velocity [23], 
i.e., 

te^S+l /2 = 2 1 + Si 8 n { U i+ 1 / 2 ) Qi + [ l ~ ^"(^+ 1 / 2 )] U i + 1/2 

= (i Qi + Q i+ 1 ) * U i + i/2 ~ \ U i + \/i\ ' fe+i - } ( 63 > 

where (i+1/2) refers to the east face of the pressure correction control volume. The above 
conventional flux estimation is only first-order accurate. The superscript (*) has been dropped in 
the above for convenience. It should be noted that the above estimation of the mass flux is identical 
to that by the AUSMD scheme (Wada & Liou [18]) except, again, the estimation of the interface 
velocity is different. 

In the present work, we also investigate a second-order estimation of the mass flux 

= <sU)Z ft + M r ?-rn) ■ ( b -m + e,-i /2 } ■ (ft - ft-,) 

+ fv^r+l/ 2 ) ' (- fc , + 3/2 + 2,+3/ 2 } ' (<?/ + 2 - ft + l) C 64 ) 

where b { _ X j 2 = U i _ l / 2 , etc. and the ratios r*_j^ 2 are the same as in Eq. (47). 

It should be noted that the remaining terms in the pressure correction equation, i.e., the ones 
involving velocity and density corrections, only contribute to the stability, not accuracy, of the 
overall algorithm. These terms vanish when overall convergence is achieved since we obtain a 
continuity-satisfying velocity field at convergence (for which pressure correction is zero 
everywhere, up to machine accuracy). 

5.3 Additional Issues Due to the Staggered Grid Layout 

An important issue in the CVS and AUSM type schemes is the estimation of interfacial 
convective velocities and pressures. Due to the staggered grid arrangement conventionally used in 
pressure— based algorithms, additional issues have to be addressed as follows. 

5.3.1 CVS 


have: 


As mentioned earlier, in the CVS, a straightforward two-point averaging is used. Thus, we 


^i+l/2 ~ ^i+l/2 J 2 (^*V + 
bj+ 1/2 = ^ij+ 1/2 = 2 ( ^'7 + ^V+t) 


(61a) 

(61b) 


It should be noted that the above interpolations result from the staggered nature of the grids 
employed for the velocity components. Due to the staggered location of the scalar variables such 
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as pressure and density, these variables are readily available at some of the control volume interfaces 
(east and west faces for u-control volumes, and north and south faces for v-control volumes). 
Wherever density is required at a point away from the scalar nodes, it is upwinded based on the 
direction of the local normal velocity component. The pressure values at the u- and v-velocity nodes 
are obtained by the pressure splitting given by Eq. (24). The local Mach numbers required for this 
splitting along the £- and ^-directions are obtained from the local normal velocity and the local 
speed of sound: 

^ = 0.5(v w+ V, v+1 )/„ 




On a staggered grid, for a u-control volume, for example, along the £- direction, the values 
of pressure at the interfaces (i+1/2 and i— 1/2) is known from the existing pressure values at the 
interfaces (due to the staggered grid, shown in Fig. 10). Pressure itself is computed from the pressure 
correction equation in which the mass fluxes are computed by upwinding density based on the 
interface convective velocity. Thus, there is no need for splitting pressure in order to get the interface 
pressures for either u— control volumes along the direction and v— control volumes along the 
^-direction. In the remaining direction for each of the control volumes, however, one can employ 
the splitting formula to obtain interface pressures. 

5.3.2 AUSM 

The AUSM type schemes use a splitting formula based on the local Mach number to compute 
the interface velocities as well as interface pressures. For a staggered grid layout, such as the one 
employed in the present algorithm, the velocity components and scalar variables (pressure and 
density) are not located at the same node. This brings up another issue (for the AUSMD scheme), 
namely, the estimation of Mach number at the locations corresponding of u, v and p, which are half 
a cell length apart from each other. For example, along the ^-direction, for the u-control volumes, 
p is located on the control volume interfaces and for the p control volumes, u is located on the 
interfaces. Thus, we need to estimate local Mach numbers along the the £— and ^—directions at the 
nodes as well as the interfaces of u, v andp control volumes. This requires, for example, the values 
of p and q (to compute the speed of sound for the Mach number computation for u-splitting) at the 
u— location and values of u (the local convection speed of for the Mach number computation for 
/^-splitting) at the p-location. Presently, this is handled via an iterative process as follows. 
Considering the ^-direction, for example 

Step 1: the Mach number at the location of p is estimated initially by Eq. (62). 

Step 2: this Mach number is used to split p and these split values of p are used to estimate p at the 
locations of u. The density values at the locations of u are obtained by upwinding based on 
the local U (normal velocity along the ^-direction). 

Step 3: Based on these p and g, the Mach number at the location of u is obtained. 

Step 4: This Mach number is used to split U and these split values of U are used to estimate an 
averaged value of U at the p locations. 

Step 5: Steps 1 through 4 are repeated a few times (typically five) until convergence is achieved. 
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6. Results of Two-Dimensional Computations 


In order to investigate the performance of the CVS and AUSM type schemes for 
two-dimensional compressible flows involving shocks, a supersonic flow over a wedge is chosen 
as a test case. Both first- and second-order CVS and AUSM schemes (minmod limiter is used for 
the second-order fluxes) for the momentum equations as well as the mass flux estimation in the 
pressure correction equation are investigated. Results obtained using the Roe scheme (implemented 
in a density-based algorithm) are also presented for comparison with the CVS and AUSM schemes. 

This problem consists of an oblique shock, generated by a supersonic flow over a wedge, 
and its subsequent reflections by a solid flat plate underneath the wedge and the wedge surface itself 
and has been investigated by Wang & Widhopf [24] among others. A schematic of the flow is 
depicted in Fig. 1 1 . The inlet Mach number is 2.9 and the wedge angle is 10.94°. Two grid systems 
are used for the computations, namely, those consisting of 101x21 and 201x41 uniformly 
distributed nodes. The location of the leading edge of the wedge is at the discontinuity in the slope 
of the top boundary of the grid layout. 

The upstream boundary condition specifies the incoming flow at the given Mach number 
whereas a zero-order extrapolation is used for the downstream boundary condition (at the exit). The 
entire bottom boundary and the wedge part of the top boundary are reflecting surfaces and thus the 
normal velocity components there are specified as zero. 

The results are presented in Figs. 12-14 in the form of thirty pressure contours with equal 
increments between the minimum and maximum pressure values. For all the cases using the CVS, 
AUSM and Roe schemes, on both the grids, the correct pressure jump and shock angles are 
predicted. However, using the first-order schemes on the 101 X 21 grid the shock is excessively 
smeared, as seen in Fig. 12. Even with the refined grid (201 x41 nodes), the first-order flux 
estimation does not yield a grid-independent solution. The accuracy improves when the momentum 
fluxes and the mass fluxes in the pressure correction equation are estimated using the second-order 
CVS along with the minmod limiter, as seen in Fig. 13. On the refined grid (201 x41 nodes), for 
example, a crisp shock structure can be observed (Fig. 13). It can be observed from the results that 
both the CVS and AUSM schemes implemented in the pressure-based solver yield accuracy 
comparable to the Roe scheme (Fig. 14). 

7. Concluding Remarks 

The separate treatment of convective and pressure fluxes is a key feature of all 
pressure-based algorithms for multi-dimensional fluid flows. Some recently developed schemes 
based on separate treatment of convective and pressure fluxes — such as the controlled variation 
scheme (CVS) the AUSM type schemes — are thus very naturally amenable for application in these 
algorithms. The approach of treating the convective and pressure fluxes in the Euler and 
Navier— Stokes equations as two distinct, though coupled, entities, appears to be very promising, as 
demonstrated by the results in the present and previous works [7,18]. The upwinding of the 
convective flux and the splitting of the pressure fluxes (based on local Mach number) achieve the 
proper propagation of signals in the system, yielding high resolution in the solution profiles with no 
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spurious oscillations. Such an approach can also be effective in the presence of source terms as 
demonstrated by the results of the longitudinal combustion instability problem. Overall, both the 
CVS and the AUSM schemes yield accuracy comparable to the Roe scheme. 

Both the CVS and the AUSM scheme yield accurate results for two-dimensional 
compressible flows, using a pressure-based algorithm. It has been demonstrated that, with these 
schemes, the pressure-based algorithms can indeed be very robust and accurate for compressible 
flows involving shocks, in addition to their well-established robustness for incompressible flows. 
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i i+1/2 i+1 


Fig. 1 . Schematic of the contributions from split pressures at an interface. 
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(b) Solution with <5 = 0.25 


Fig. 2. 


Total energy profiles for the shock tube problem using Formulation 1 of the CVS ip term 
only in the pressure flux) with two values of <3; minmod limiter is used. 
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(b) Solution with <5 = 0.25 
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Total energy profiles for the shock tube problem using Formulation 2 of the CVS (p and 
pu terms in the pressure flux) with two values of <3; minmod limiter is used. 
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(b) Formulation 2 (p and pu terms in the pressure flux) 

Fig. 4. Total energy profiles for the shock tube problem using the second-order AUSM scheme 
with the two formulations for the pressure flux; minmod limiter is used. 
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(a) Solution with <5 = 0.0 
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(b) Solution with <5 = 0.5 

Total energy profiles for the shock tube problem using Formulation 1 of the CVS with 
the pressure splitting based on interface Mach number, i.e., Eq. (27) (minmod limiter). 




Fig. 6. Total energy profile using the second— order AUSM scheme for the convective flux 

(minmod limiter) with pressure splitting based on interface Mach number, i.e., Eq. (27). 



Fig. 7. Total energy profiles for the shock tube problem using the second— order Roe scheme, 
minmod limiter is used. 
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(iii) Roe scheme (iii) Roe scheme 

(a) Pressure (b) Temperature 


Fig. 8. Ten pressure and temperature modeshapes for the combustion instability problem 
using second-order CVS, AUSM and Roe schemes (with minmod limiter). 
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PRESSURE PRESSURE PRESSURE 




(a) Illustration of the staggered locations of u, v and p, and the nomenclature for a typical 

H-control volume. 

(*.»/+ U+i 



(b) Notation of the staggered grid in curvilinear coordinates. 

Fig. 10. Location of variables u, v and p on a staggered grid for the pressure-based algorithm. 
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(a) First-order CVS for all the equations on the 101 X 21 grid. 



(b) First-order AUSM for all the equations on the 101 X21 grid. 



(c) First-order CVS for all the equations on the 201 X41 grid. 



(d) First-order AUSM for all the equations on the 201 x41 grid. 


Fig. 12. Pressure contours for a supersonic flow (inlet Mach number = 2.9) over a wedge (angle 
10.94°) on the 101 X21 and 201 x41 grids using first-order CVS and AUSM schemes. 
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(a) Second-order CVS for all the equations on the 101 x 21 grid. 



(b) Second-order AUSM for all the equations on the 101 x21 grid. 



(c) Second-order CVS for all the equations on the 201 X41 grid. 


(d) Second-order AUSM for all the equations on the 201 X 41 grid. 



Fig. 13. Pressure contours for a supersonic flow (inlet Mach number = 2.9) over a wedge (angle 
10.94°) on the 101 x21 and 201 X41 grids using second-order CVS and AUSM schemes. 
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Fig. 14. Pressure contours for a supersonic flow (inlet Mach number = 2.9) over a wedge (angle 
10.94°) on the 101 x21 and 201 X41 grids using the second-order Roe scheme. 


37 



REPORT DOCUMENTATION PAGE 


Form Approved 
OMB No. 0704-0188 


Public reporting burden for this collection of information Is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 
gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 
collection of information. Including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports. 1215 Jefferson 


Davis Highway. Suite 1204, 


A 22202*4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-01 88). Washington. DC 20503. 


1. AGENCY USE ONLY ( Leave bianK) 


4. TITLE AND SUBTITLE 


2. REPORT DATE 

March 1995 


3. REPORT TYPE AND DATES COVERED 

Technical Memorandum 


Investigation of Convection and Pressure Treatment with Splitting Techniques 


6. AUTHOR(S) 


Siddharth Thakur, Wei Shyy, and Meng-Sing Liou 


5. FUNDING NUMBERS 


WU-505-90-5K 


7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135-3191 


8. PERFORMING ORGANIZATION 
REPORT NUMBER 


E-9483 


9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESSES) 

National Aeronautics and Space Administration 
Washington, D.C. 20546-0001 


10. SPONSORING/MONITORING 
AGENCY REPORT NUMBER 

NASA TM- 106868 
ICOMP-95-3 


11. SUPPLEMENTARY NOTES 

Siddharth Thakur, University of Florida, Department of Aerospace Engineering, Mechanics and Engineering Science, Gainesville, Florida 32611; Wei 
Shyy, Institute for Computational Mechanics in Propulsion, Lewis Research Center, Cleveland, Ohio and University of Florida, Department of Aerospace 
Engineering, Mechanics and Engineering Science, Gainesville, Florida 32611 (work funded under NASA Cooperative Agreement NCC3-370); and 
Meng-Sing Liou, NASA Lewis Research Center. ICQMP Program Director, Louis A. Povinelli, organization code 2600, (216) 433-5818. 


12a. DISTRIBUTION/ AVAILABILITY STATEMENT 12b. DISTRIBUTION CODE 

Unclassified - Unlimited 
Subject Categories 34 and 64 

This publication is available from the NASA Center for Aerospace Information, (301) 621-0390. 


13. ABSTRACT (Maximum 200 words) 

Treatment of convective and pressure fluxes in the Euler and Navier-Stokes equations using splitting formulas for 
convective velocity and pressure is investigated. Two schemes — Controlled Variation Scheme (CVS) and Advection 
Upstream Splitting Method (AUSM) — are explored for their accuracy in resolving sharp gradients in flows involving 
moving or reflecting shock waves as well as a one-dimensional combusting flow with a strong heat release source term. 
For two-dimensional compressible flow computations, these two schemes are implemented in one of the pressure-based 
algorithms, whose very basis is the separate treatment of convective and pressure fluxes. For the convective fluxes in the 
momentum equations as well as the estimation of mass fluxes in the pressure correction equation (which is derived from 
the momentum and continuity equations) of the present algorithm, both first- and second order (with minmod limiter) flux 
estimations are employed. Some issues resulting from the conventional use in pressure-based methods of a staggered grid, 
for the location of velocity components and pressure, are also addressed. Using the second-order fluxes, both CVS and 
AUSM type schemes exhibit sharp resolution. Overall, the combination of upwinding and splitting for the convective and 
pressure fluxes separately exhibits robust performance for a variety of flows and is particularly amenable for adoption in 
pressure-based methods. 


14. SUBJECT TERMS 


Convective fluxes; Pressure fluxes; Flux splitting 


17. SECURITY CLASSIFICATION 18. SECURITY CLASSIFICATION 19. SECURITY CLASSIFICATION 
OF REPORT OF THIS PAGE OF ABSTRACT 

Unclassified Unclassified Unclassified 


15. NUMBER OF PAGES 

39 


16. PRICE CODE 

A03 


20. LIMITATION OF ABSTRACT 


NSN 7540-01-280-5500 


Standard Form 298 (Rev. 2-89) 

Prescribed by ANSI Std. Z39-18 
298-102 
























